The wall shear rate distribution for flow in random sphere packings 
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The wall shear rate distribution P(y) is investigated for pressure-driven Stokes flow through 
random arrangements of spheres at packing fractions 0.1 < (f> < 0.64. For dense packings, P{^) is 
monotonic and approximately exponential. As <f> — » 0.1, P(j) picks up additional structure which 
corresponds to the flow around isolated spheres, for which an exact result can be obtained. A simple 
expression for the mean wall shear rate is presented, based on a force-balance argument. 

PACS numbers: 47.56.+r, 47.15.G- 



The wall shear rate 7 is the rate at which the tangen- 
tial velocity of a fluid vanishes on approaching a wall. 
It determines the hydrodynamic forces acting on a par- 
ticle adjacent to the wall, and is therefore a key quan- 
tity governing deposition, retention, and detachment [l| . 
For example, 7 is an important factor which determines 
whether colloidal particles can become attached to a sur- 
face by specific ligand binding 0]. Experimentally, such 
processes are often examined using flow cells with well- 
controlled hydrodynamics for which the wall shear rate 
is known. For flow through a porous material though, 
one can expect a distribution of wall shear rates P("f)- 
This is illustrated in Fig. Q] A crucial issue for realis- 
tic situations, such as deep bed filtration 0] or partic- 
ulate soil detergency in fabric cleaning 0], is therefore 
to characterise the wall shear rate distribution ^(7) for 
flow in more complex pore spaces. This is also a prob- 
lem of generic interest in the growing field of statistical 
microhydrodynamics. Previously, only P{^) for flow in 
two-dimensional channels with (fractally) rough walls has 
been investigated [jj. 

In this Letter, P(7) and its relation to the mean fluid 
velocity U m is investigated for pressure-driven Stokes flow 
in random sphere packings at packing fractions in the 
range 0.1 < (j> < 0.64. The relationship between P(~f) 
and U m is of paramount importance for applications since 
it is very difficult to access P(j) experimentally (it either 
has to be done by detailed resolution of the flow field, 
or indirectly by looking at the behaviour of particulate 
tracers), but determination of U m is much easier. 

We generated sphere packings in periodic simulation 
boxes for packing fractions in the range 0.1 < <fi < 0.64 
by a Monte-Carlo algorithm [f|. The highest packing 
fraction corresponds to the usual random close-packed 
limit. Whilst the lower packing fractions are mechan- 
ically unstable, they provide a useful interpolation be- 
tween isolated spheres and packed beds. We also gen- 
erated a slightly looser packing of touching spheres at 
(j) 1=3 0.622 by a sequential deposition algorithm [?]. This 
latter geometry is not periodic in the z-direction (the 
deposition direction), but we have found that the bulk 
properties can be determined by extrapolation. 



For the flow-field calculations we use a standard lat- 
tice Boltzmann (LB) methodology which is now well- 
developed for this class of problems 0,0,13, 11, 12, 13, 
[Hl,[l5j]. A s already mentioned, we solve the Stokes equa- 
tions and thus operate at strictly zero Reynolds number. 
The spheres are held stationary and flow is generated in 
the pore space by applying a uniform body force corre- 
sponding to a mean pressure gradient (Vp) m in the x-, y- 
or z-directions. The hydrodynamic forces exerted on wall 
lattice nodes are easily found in LB. For each wall node 
one can trivially determine the tangential force compo- 
nent since the corresponding sphere centre is known. The 
local wall shear rate is then given by the magnitude of 
the tangential component divided by the viscosity. In 
this way we obtain a large set of wall shear rates from 
which we reconstruct P{^) [16]. We also measure the 
mean volumetric (superficial) fluid velocity U m . 

We first discuss our results for the permeability, k, 
since this underpins our analysis of P{"i). It is defined via 
Darcy's law, U m — (k/ri)(S7p) m , where 77 is the viscos- 
ity. Our results, expressed in dimensionless terms using 
k/a 2 , are shown as a function of packing fraction in Ta- 
ble U and Fig. [2] Generally speaking, the permeability 
falls dramatically with increasing packing fraction. For 
4> < 0.5 our results are in excellent agreement with pre- 
vious work by Ladd [ItJ and van der Hoef et al. [IH ■ For 
(j) > 0.6 our results are ~ 10% higher than the accurate 
results obtained recently by van der Hoef et al. [l3j], al- 
though we are in agreement with some previous studies 
[2, HH . This may reflect subtle differences in the way the 
sphere packings are constructed. The sequential deposi- 
tion packing at <fi « 0.622 fits nicely into the series. In 
this case the permeability is in principle different parallel 
to and perpendicular to the deposition direction. We find 
though that the difference is certainly less than 10%, in 
agreement with Coelho et al. 0]. 

An oft-used correlation is the Kozeny-Carman relation, 



k = (1 - & 3 /c s 2 



(1) 



where s = 64>/o~ is the specific surface area of spheres 
(in any arrangement ) and the numerical factor cq ~ 4- 
5 0, [12, 0, Ell ■ We find this does indeed capture the 
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FIG. 1: The flow in the controlled geometry of a flow cell gives 
rise to a uniform wall shear rate (left), whereas the flow in a 
porous material gives rise to a distribution of wall shear rates 
(right). It is the wall shear rate 7 that governs the deposition 
and detachment of particles (inset). 

behaviour of the permeability quite well for intermedi- 
ate to high packing fractions (Table [!}. Interestingly, 
for tj> > 0.2 we noticed our data can be accurately fit by 
log(A;/cr 2 ) = A+B(f> with A = —1.04(6) and B — -9.6(1), 
reminiscent of what has been found for fibrous beds [l0( • 
Now we turn to the mean wall shear rate, defined via 
7m = /q ^7 7-f(7) • F° r Stokes flow, j m is strictly pro- 
portional to U m , so that <J~f m /U m is a convenient way to 
express the mean wall shear rate in dimensionless terms, 
shown in Table Hand Fig. O We see that o-j m /U m grows 
dramatically with packing fraction, similar to the inverse 
of k/a 2 . 

This behaviour can be understood by the following 
force-balance argument. The force per unit volume act- 
ing on the fluid due to the mean pressure gradient is 
(1 — 4>)(Vp) m . In steady state this must balance the in- 
tegrated wall stress, thus the mean wall stress is exactly 
(1 — </>)(Vp) m / s where s is the specific surface area. If 
we now approximate the mean wall stress by ?77 m , use 
Darcy's law to replace (Vp) m by U m , and substitute 
s = 6(/>/cr, we get 

7m = a(l - <f>)aU m /(6<f>k). (2) 

We have introduced a prefactor a to capture the approx- 
imate nature of this expression. From our data we find 
that a ~ 0.6-0.8 is very insensitive to packing fraction, 
as shown in Table [J (we can rationalise this value of a by 
arguing that, on average, 2/3 of the wall stress lies in the 
wall tangent plane). Eq. @ explains the approximate 
inverse relationship between aj m /U m and k/a 2 . Inci- 
dentally, in a parallel-sided capillary of arbitrary cross 
section, the flow is laminar and parallel to the walls. In 
this case the mean wall stress is exactly rjj m and Eq. (|2|) 
is exact with a = 1. Our LB methodology is constructed 
to retain this exact result, provided the capillary axis is 
aligned with a grid axis. 

Finally we turn to the wall shear rate distribution, 
which we report in terms of x — "f/"/ m and f(x) defined 
such that P( 7 ) = (l/7m) /(7/7m)- At packing fractions 
4> > 0.6, f(x) is monotonic and quite well approximated 
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TABLE I: Dimensionless permeability k/a 2 and mean wall 
shear rate a^ m /U m as a function of packing fraction 0: Co is 
the Kozeny-Carman factor in Eq. (fl]) and a is the prefactor 
in the force-balance expression in Eq. For the sequential 
deposition sample (</> ~ 0.622), results are given parallel and 
perpendicular to the deposition direction (z). A figure in 
brackets is an estimate of the error in the final digit [15|] . 

by an exponential (Fig. [31 upper plot). It is interesting to 
note that a similar exponential distribution is found for 
the local flow speeds although in this case a peak at zero 
is to be expected given the large volume of pore space 
immediately adjacent to the sphere surfaces [ll|, EH- We 
will return to the small x behaviour of f(x) in a moment. 

As the packing fraction is reduced, a hump appears in 
f(x) at around x = 0.5-0.6 (Fig.[3J lower plot). This fea- 
ture seems to be associated with the transition from chan- 
nel flow at high packing fractions towards flow around 
individual spheres at lower packing fractions. This in- 
terpretation is supported by the exact result which can 
be obtained for P(7) from the Stokes solution for flow 
around a sphere, as we now discuss. 

A remarkable feature of Stokes flow around a sphere is 
that the wall stress has the same vectorial value 3r/U/cr 
at all points on the sphere surface, where U is the flow 
velocity at infinity [20j. If we project this into the 
wall tangent plane, we obtain the local wall shear rate 
7 = (3U m sin#)/cr, where 8 is the angle between the wall 
normal and the direction of the flow field at infinity, and 
U m = |U|. The mean wall shear rate is then given by 
ajm/Um = J* (3/2)sm 2 0d9 = 3tt/4 w 2.356. It fol- 
lows that x — 7/7m = (4/7t) sin 9, and from f(x)dx = 
(1/2) sin^ d9 (i. e. the area measure fig). 

f(x) = Q<x<A/tt. (3) 

V ; V(4/tt) 2 -x 2 - - I ^ > 

This is the desired exact result for the wall shear rate 
distribution for Stokes flow around an isolated sphere, 
shown as the dotted line in the lower plot of Fig. [3] It 
diverges as x — > 4/7T w 1.273, corresponding to 9 — > 7r/2 
where the wall shear rate is maximal. This behaviour 
is, we believe, responsible for the hump that appears 
in f(x) at low packing fractions. The fact that there 
is still a significant difference between Eq. ^ and f(x) 
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FIG. 2: Dimensionless permeability and mean wall shear rate 
as a function of packing fraction, from Table [I] The solid line 
is Eq. (31) from van der Hoef et al. [l3T ] which is claimed to 
be accurate to within 3%. The dashed line for the mean wall 
shear rate data is a guide to the eye. Error bars are smaller 
than the symbols. 



for = 0.1 should not be too surprising given the long 
range nature of hydrodynamic interactions. We see this 
also in k and 7 m which are, respectively, a factor « 2.76 
smaller and a factor « 1.9 higher, than the correspond- 
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ing isolated sphere limits (i. e. k/a 2 — 1/(180) 
and crj m /U, n — 37r/4 derived above). In fact the perme- 
ability data from Ladd suggests that the isolated sphere 
result is approached only very slowly as cj> — > [TtJ • 

Now we return to the small x behaviour of f(x). 
Clearly, for any sphere, the local wall shear rate has 
to vanish at least at one point on the sphere surface — 
this is a consequence of the so-called 'hairy ball theo- 
rem' [22|. Thus it is not at first sight surprising that 
f(x) goes to a plateau as x — > (Fig. [31 lower plot). 
However, Eq. ([3]) has the property that f(x) ~ i as 
x — > arising from the stagnation points at 9 = (0, w). 
This behaviour might be expected to be generic for low 
packing fractions where stagnation points are common. 
In contrast, for dense sphere packings the flow is more 
channel-like and stagnation points are rare. In this case 
the wall shear rate vanishes, inter alia, at all the con- 
tact points between spheres. Analysis of pressure-driven 
flow in the vicinity of a contact point using the Reynolds 
lubrication equation Q suggests f(x) ~ x 5 for x — > 
where 5 = (4- > /10)/(v'10 - 2) w 0.72. It is therefore 
rather surprising that, independent of packing fraction, 
a plateau rather than a power law is observed for f(x) as 
x -> 0. 

One possible reason for this is that long-range flow field 
inhomogeneities (on length scales > a) wash out the ex- 
pected behaviour and replace the power law by a plateau. 
We investigated this possibility by constructing an indi- 
vidual f(x) for each sphere, then averaging over all the 
spheres in a sample. This should remove the effects of 
long-range flow field inhomogeneities. We find though 




FIG. 3: The upper plot shows the wall shear rate distributions 
for all data sets with (j> > 0.6. The dashed line is f(x) — e~ x . 
The lower plot shows the same for the six periodic packings 
with 0.1 < (j> < 0.64; the curves are displaced for clarity. The 
dotted line is the exact result in Eq. (|3} for Stokes flow around 
an isolated sphere. 



there is little change in f(x); the hump at low (f> becomes 
somewhat more pronounced but the plateau remains in 
all cases. At the same time we also examined the hydro- 
dynamic forces acting on individual spheres. We found 
that these have a relatively narrow distribution (approx- 
imately Gaussian, with a standard deviation 20-30% of 
the mean) indicating that the flow field on length scales 
> a is rather homogeneous. We conclude that long-range 
flow field inhomogeneities are unlikely to be important. 
Instead, the implication is that the shape of f(x), and in 
particular the plateau at x — ► 0, is mostly controlled by 
the local pore geometry. The important message seems 
to be that using highly idealised situations, such as the 
Stokes solution for flow around an isolated sphere or lu- 
brication theory in the vicinity of a contact point, may 
give qualitatively misleading results when it comes to in- 
ferring the overall statistical properties. 

To summarise, for applications Eq. |2]) provides the 
key link between the mean wall shear rate j m and the 
mean fluid velocity U m . If necessary the Darcy per- 
meability can be estimated from the Kozeny-Carman 
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relation in Eq. ([T]). Knowledge of j m is then suffi- 
cient to determine the whole wall shear rate distribu- 
tion, if the latter is assumed to be exponential, i. e. 
P(j) 1=3 (l/7 m ) exp(— 7/7 m ). More generally, our investi- 
gation demonstrates how direct numerical calculation of 
the statistical properties of microhydrodynamic flows can 
complement exact solutions for simplified geometries, to 
gain new insights. 
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